NPS71 -88-001 

NAVAL 



(REVISED) 

POSTGRADUATE SCHOOL 

Monterey, California 




A PROGRAM TO COMPUTE MAGNETIC ANOMALY 
DETECTION PROBABILITIES 



R. N, FORREST 



OCTOBER 1988 

Approved for public release; distribution is unlimited. 

Prepared for: 

Naval Postgraduate School 
Monterey, CA 93943-5000 



FEDDOCS 
D 208.14/2 
NPS-71-88-001/ 
988 




I 



'■i 1. 



' f> \ 



NAVAL POSTGRADUATE SCHOOL 
Monterey, California 



Rear Admiral R. C. Austin Harrison Shull 

Superintendent Provost 

Reproduction of all or part of this report is authorized. 

This report was prepared by: 



Unclassified 



Security Classification of this page 



REPORT DOCUMENTATION PAGE 



la Report Security Classification Unclassified 



lb Restrictive Markings 



2a Security Classification Authority* 



2b Declassification/Downgrading Schedule 



3 Distribution Availability of Report 

Approved for public release; distribution is unlimited. 



4 Performing Organization Report Number(s) 

NPS71-88-001 (Revised) 



5 Monitoring Organization Report Number(s) 



6a Name of Performing Organization 

Naval Postgraduate School 



6b Office Symbol 
(If Applicable) 1 1 



7a Name of Monitoring Organization 

Naval Postgraduate School 



6c Address (city, state, and ZIP code) 

Monterey, CA 93943-5000 



7 b Address (city, state, and ZIP code) 

Monterey, CA 93943-5000 



8 a Name of Funding/Sponsoring Organization 



8b Office Symbol 
(If Applicable) 



9 Procurement Instrument Identification Number 



DIRECT FUNDING 



8c Address (city, stale, and ZIP code) 



1 0 Source of Funding Numbers 

Program Element Number | Project No [ Task No | Woric Unit Acccssicffi 



1 1 Title (Include Security Classification) 

A PROGRAM TO COMPUTE MAGNETIC ANOMALY DETECTION PROBABILITIES 



12 Personal Author(s) R. N. FoirCSt 



13a Type of Report 


13b Time Covered 


14 Date of Report (year, month, day) 


15 Page Count 


Technical 


From To 


October 1988 


70 



1 6 Supplementary Notation 



17 Cosati Codes 


Field 


Group 


Subgroup 









18 Subject Terms (continue on reverse if necessary and identify by block number) 

Magnetic Anomaly Detection, MAD, Lateral Range Function, Detection Models, 
Magnetic Dipole 



19 Abstract (continue on reverse if necessary and identify by block number) 

The report contains user instructions, a listing and documentation for a microcomputer BASIC program that 
can be used to compute an estimate of the probability that a magnetic anomaly detection (MAD) system such as the 
AN/ASQ-81 will detect a submarine during an encounter. 



20 Distribution/Availability of Abstract 

|x| unclassified/unlimited QJ same as report QJ DTICi 



21 Abstract Security Classification 

Unclassified 



22a Name of Responsible Individual 

R. N. Forrest 



22b 



Telephone (Include Area code) 

(408) 646-2653 



22c Office Symbol 

55Fo 



DD FORM 1473. 84 MAR 



83 APR edition may be used until exhausted 
All other editions are obsolete 



security classification of this page 

Unclassified 



1 



The program in this report is presented without representation or 
warranty of any kind. 



Table of Contents 



I. Introduction 1 

II. Program User Instructions 3 

III. Encounter Model Limitations 6 

IV. Parameter Values 9 

Appendix 1. The Detection Statistics 12 

Appendix 2. Program Probability Calculations 18 

Appendix 3. The Magnetic Signal 20 

Appendix 4. The Anderson Formulation 24 

Appendix 5. The Encounter Equations of Motion 30 

Appendix 6. The Submarine Magnetic Dipole 3 3 

Appendix 7. The Earth Magnetic Field and Dip Angle- 36 

Appendix 8. An Alternative Encounter Model 41 

Appendix 9. An Example of the Program Output 4 3 

Appendix 10. A Program Listing 54 

References 64 



I. Introduction 

This report contains user instructions, a listing and 
documentation for a microcomputer BASIC program that can be used 
to compute an estimate of the probability that a magnetometer 
based detection system such as the AN/ASQ-81 will detect a 
submarine magnetic dipole field during an encounter. 

The program generates detection probabilities based on two 
encounter models. In the first encounter model, the detection 
system uses a crosscorrelation detector. In the second 
encounter model, the detection system uses a square law detector. 
Relative to operationally realizable values, probabilities based 
on the first model represent upper bounds and those based on the 
second model represent lower bounds. For both encounter models, 
the signal is proportional to the square of the magnitude of the 
projection of the dipole field on the earth magnetic field and 
the magnetic noise does not change with changes in the position 
of the magnetometer. Also, for both encounter models, an 
encounter is a straight line encounter with constant vertical 
separation. 

The encounter models can be interpreted as models of a 
magnetic anomaly detection (MAD) system on an aircraft that is 
moving with constant course, speed and altitude in an encounter 
with a submarine moving with constant course, speed and depth. 

Or, they can be interpreted as models of a stationary magnetic 
anomaly detection system in an encounter with a submarine moving 
with constant course, speed and depth. 
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The program parameters include: encounter latitude and 

longitude, submarine induced magnetic moments, submarine 
permanent magnetic moments, submarine course, speed and altitude, 
magnetometer course, speed and depth, encounter lateral range 
(the horizontal range at the closest point of approach in a 
straight line encounter) and false alarm rate. (In the program, 
a false alarm is the event that the detection system classifies 
noise as a dipole signal.) 
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II. Program User Instructions 

As listed in Appendix 10, the program can be run under a 
BASIC language that is compatible with IBM PC BASIC. If the 
listing is used to enter the program through the keyboard, then 
the program should be saved with the name MAD. BAS or the value of 
N$ on line 40 should be changed to the file name under which it 
is saved. 

The program contains user instructions in the form of query 
and parameter limitation messages. As an example of the former, 
after starting the program under BASIC, the following message 
should appear: 

Magnetic Anomaly Detection (MAD) Lateral Range Function 

generate or print a program data file (g/p)? 

By entering p, data can be printed from a program data file 
that was generated by the program. By entering g, a program 
data file can be generated for a set of user specified 
conditions. With either response, a sequence of additional 
queries is displayed. These queries require either an indicated 
response or a parameter value as the input. If the initial 
response is g, the sequence includes queries whose responses 
determine whether or not an auxiliary data file will be generated 
that can be used for future input of magnetic, processing or 
kinematic data. In particular, the first query in the sequence 
gives the option of using a combined magnetic, processing and 
kinematic data file. If the response indicates that it should 
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be used and the file is available, the parameter values that 
remain to be entered in order to generate a program data file are 
the following: the false alarm rate, the magnetic noise, the 
maximum encounter lateral range and the lateral range step. The 
combined file should be used only if the effect of varying just 
one or more of these parameters is desired. If the response 
indicates that the file will not be used, queries concerning 
magnetic, kinematic and processing parameter values are 
displayed. 

After all of the program parameter values have been 
entered, a query is displayed giving the option of generating the 
combined file. Then a query giving the option of generating a 
program data file, a query giving the option of printing the 
encounter parameter values and a query giving the option of 
printing lateral range function values are displayed. The 
lateral range function values are the encounter detection 
probabilities indexed by lateral range. The parameters maximum 
lateral range and lateral range step determine the index lateral 
ranges of the encounters for which probabilities of detection are 
computed . 

The program generates magnetic signal values that correspond 
to points in time during an encounter. Following the lateral 
range function query, a query is displayed that gives the option 
of printing the magnetic signal values for an encounter. If the 
option is exercised, the option is repeated. When the option is 
not exercised, a query is displayed that gives the option of 
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generating or printing a new program data file. If this option 
is not exercised, the program ends. 

Some suggested guides for determining parameter values can 
be found in Section IV of this report. 
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III. Encounter Model Limitations 

In the two encounter models, the signal is the unfiltered 
signal that is generated by a magnetic dipole that moves relative 
to a magnetometer. Depending on the input filtering, describing 
a ship magnetic anomaly as a dipole anomaly should not be a 
significant limitation for an encounter slant range at the 
closest point of approach (CPA) that is greater than one hull 
length. The detection decision is based on samples from a single 
time interval (window) that is centered on the CPA. The length 
of the sampling interval and the sampling rate are parameter 
values that are inputs to the program. In terms of signal-to- 
noise ratio, there is an optimum sampling interval length 
(integration time) and sample rate. Although dipole signal 
energy is not symmetrically distributed about the CPA time, for a 
given sampling interval, the difference between the signal energy 
for the optimum interval location and the CPA centered location 
should not be significant in most cases. 

The encounter model magnetometer noise samples are 
determined by a gaussian random process. They are values of 
independent identically distributed normal random variables. The 
standard deviation of these random variables is referred to as 
the magnetic noise and the variance is the noise in the sense of 
the signal-to-noise ratio. The degree of correspondence between 
this process and operational noise depends on the nature of the 
dominant operational noise sources and on the magnetometer input 
filtering (noise whitening) . 
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In the encounter models, the intervals are adjacent but not 
overlapping and a detection statistic corresponds to each sample 
interval and its value is determined by the sample values. If 
the value of the detection statistic for an interval equals or 
exceeds a threshold value, a detection is indicated. The 
threshold value is determined by the false alarm probability 
which in turn is determined by the false alarm rate and the 
sampling interval length. 

A detector that used a moving sample interval that was 
generated by replacing the oldest sample by the newest one would 
correspond more closely to the detector in an operational 
detection system. In a model of the detector, since the sample 
windows overlap, the detection statistics would represent a 
sequence of random variables that were correlated over an 
interval equal to the width of the sample interval. Because of 
this dependence, it seems unlikely that the results that would be 
obtained with an encounter model based on an overlapping interval 
detector would differ significantly from those obtained with an 
encounter model based on a nonoverlapping interval detector. 

For encounter lengths of the order of a few nautical miles 
or less, the straight line encounter condition should not be a 
significant limitation. In particular, this should be the case 
for a fixed magnetometer since, for a submarine (or surface ship) 
target, vertical separation and course changes should be less 
likely to occur. 
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other models are available that can be used as the basis for 
computing an estimate of the probability that a magnetic anomaly 
detection system will detect a submarine during an encounter. 

For example, one that is described in Appendix 8 can be used to 
determine the slant ranges of straight line encounters for which 
the detection probability is equal to a specified value. The 
parameter values that are recpaired to do this are an average 
submarine dipole moment, a detection system capability factor and 
a noise factor. Values for these parameters can be determined 
from operational data. However, the values are specific to 
averages over a particular set of encounter conditions. An 
advantage of the two encounter models relative to this model is 
their adaptability to different magnetic, processing and 
kinematic conditions. 
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IV. Parameter Values 



The magnetic parameter value queries are generally 
explanatory with regard to the value that should be entered. 

This is also true of the kinematic parameter value queries. 
However, there is some ambiguity with respect the processing 
parameter value queries and the noise parameter value query. To 
reduce this ambiguity, a brief discussion of the common 
characteristics of the two encounter model processing parameters 
and the noise parameter is given below. This is followed by 
some guidelines for choosing these parameter values. 

In both encounter models, a decision is made at the end of 
each sampling interval. The decision is either noise energy was 
present during the interval or noise energy plus signal energy 
was present during the interval. The sample intervals are 
adjacent, equal width, nonoverlapping time intervals. The number 
of samples that are input in a sample interval is determined by 
the sampling rate and the interval length. 

The program default choice for the sampling rate is 2-MAXF 
where MAXF is a parameter that is labeled the maximum magnetic 
signal frequency. This sample rate is the Nyquist rate for an 
ideal low pass filter. However, the signal in a sample interval 
that is computed by the program represents an unfiltered dipole 
signal. This is a reasonable approximation if the signal energy 
that is associated with signal components greater than MAXF is 
relatively small. As discussed below, the noise energy in a 
sample interval should be considered to be proportional to MAXF 
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in order to be consistent with the encounter models. Ideally, a 
default choice for MAXF would make the ratio of the signal energy 
to the noise energy a maximum for the sampling interval of an 
encounter. The program default choice for MAXF is 2 • MAXVM/MINRO 
where MAXVM is a user estimated maximum encounter relative speed 
converted from knots to meters per second and MINRO is a user 
estimated maximum slant range at CPA in meters in terms of a just 
detectable target. (A precise definition of a just detectable 
target can be made in terms of a specified detection probability, 
false alarm probability and target dipole moment.) The default 
choice for MAXF is consistent with the observation in Reference 1 
that if an optimum value for MAXF is determined for a minimum 
dipole moment target, then no significant increase in MAXF is 
required in order to maintain a required detection probability if 
the encounter lateral range is decreased even though the signal 
energy spectrum is shifted to higher frequencies. 

Because of the detection statistics and the gaussian noise 
model, if there were no penalty for decision delay, a sample 
interval length for a signal should be chosen equal to the signal 
duration, since this would make the detection probability for an 
encounter a maximum. In the program, a default choice for the 
sample interval length (integration time) is 2 • MAXRO/MINVM where 
MINVM is a user estimated maximum encounter relative speed 
converted from knots to meters per second and MAXRO is a user 
estimated maximum slant range at CPA in meters for a detectable 
encounter (in terms of a specific detection probability and false 
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alarm probability) with a user specified maximum dipole moment 
target. This choice might be considered a balance between 
minimizing decision delay and maximizing detection probability. 

In the encounter models, the sample intervals are located so that 
the CPA time is at the center of a sample interval and this 
sample interval is the only one that contains signal energy. 

This characteristic is consistent with the default choice for the 
sample interval. For a given sample interval length, although in 
general the interval location is not the optimum one in terms of 
signal energy, it should be approximately so in most cases. 

The program noise parameter is SIG. It represents the 
standard deviation a associated with the magnetic noise process 
of the two encounter models. In terms of the ideal low pass 
filter implied by the encounter models, its square should be 
equal to MAXF- (SIGO) ' where SIGO is the magnetic noise process 
constant power spectral density. The program does not enforce 
this relation. Therefore, in using the program, the implied 
relation between the two input parameters: magnetic noise and 
maximum magnetic signal frequency should be kept in mind. If an 
average value of the peak-to-peak magnetic noise for an encounter 
can be estimated, for example from a magnetometer trace, then the 
value for SIG should be chosen so that the estimate is 4 to 6 
times this value. 



11 



Appendix 1. The Detection Statistics 

In this appendix, Y2> ' ' ' > sequential values 

(voltages) representing the sample values in a sample interval. 
They are the input to a magnetometer's detector. With these m 
sample values, the detector computes the value of a detection 
statistic. This value is represented by x and the detection 
decision corresponding to the sample interval is determined by 
the decision rule: If x > x*, then the input during the sample 
interval was noise plus signal, otherwise, the input was noise. 

For both encounter models, the detection probability and the 
false alarm probability are decreasing functions of x and the 
relation is one-to-one in both cases. In the program, the false 
alarm probability pf is used to determine a unique value of 
the threshold x*. This value is then used to determine a unique 
value of the detection probability p^j. 

In the program, pf is found using the relation pf = R-5t 
where R is the false alarm rate in false alarms per second 
and St is the sample interval length in seconds. This relation 
is based on the following argument: With no signal energy in a 
sample interval, y^^ = n^^, y 2 = n 2 / • • • / Ym “ ^m where n^^, n 2 , 

• ■ ■ , ni,^ are noise values (voltages) input to a magnetometer's 

detector. The n^ are values of independent normal (gaussian) 
random variables, each with mean zero and standard deviation a. 
Because of this, in the encounter models, values of x for 
different sample intervals are the values of independent random 
variables that determine two outcomes: x > x* or x < x*. 
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Therefore, in terms of these outcomes, an encounter is a sequence 
of independent Bernoulli trials. If there is no signal present, 
since the noise is determined by a stationary process, pf is the 
same for each sample interval and, in terms of these outcomes, 
the sequence is a series of repeated independent Bernoulli trials 
and therefore the expected number of trials between false alarms 
is 1/Pf. Since the time between trials is 6t, the expected 
number of seconds between false alarm is 6t/pf or the expected 
number of false alarms per second R is equal to pf/6t. 

The determination of x* depends on the encounter model 
statistic. For both encounter models, when there is a signal, 
yi = ni + si, y 2 = n 2 + y 2 , • • • / Ym = >^m + where s^, S 2 , 

• • • , are signal values (voltages) input to a magnetometer's 

detector. The models imply that the signal values = K- (Hg) i 

where K is a constant whose value is determined by the 
characteristics of the encounter magnetometer and where the 
(Hs) i are dipole magnetic signal intensities. The models also 
imply that the noise values nj^ are determined by a gaussian 
stochastic processes characterized by a standard deviation a 
and that nj^ = + nj^ where the (Hj^) are magnetic noise 

intensities and the n| are magnetometer instrument noise values 
(voltages). In the program, the magnitude of K is 1. Since, 
for both models, P (3 depends only on the ratio of signal energy 
to noise energy for a sample interval, this is a satisfactory 
choice for the program if instrument noise is assumed to be 
determined by a process independent of the magnetic noise process 
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and to be expressed in terms of an equivalent magnetic noise 
(H)l = (1/K) -nl. 

For both encounter models, the signal (the average signal 
power) S = (1/m) -E s^ and the noise (the expected value of the 
average noise power) N = a so that the signal— to— noise ratio 
is (1/m) -Z s[/o^ where the sum index i = 1, 2, • • • , m. 

The Crosscorrelation Detector Statistic: The statistic for the 

first encounter model is a crosscorrelation detector statistic 
that is defined by the sum 

X = Z Yi'Si 

where the summation index i = 1, 2, • • • , m and the sum is over 

the values corresponding to a sample interval. For the first 
encounter model, the characteristics of both the noise and the 
signal are required in order to determine encounter detection 
probabilities. In particular, the signal values for an encounter 
are in the memory of the detector prior to the encounter. For 
the encounter conditions and a specified false alarm 
probability, the statistic is optimum in the sense that the 
encounter detection probability for this statistic is at least 
equal to that for any other statistic. Because of these 
considerations, encounter probabilities based on the 
crosscorrelation statistic can be considered to represent upper 
bounds on detection performance against dipole targets for 
magnetometers of the type described by the models. 

For a sample interval without signal energy, x is the value 
of a normal random variable with a mean = 0 and a variance 
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2 2 2 , , , 

= a -S Sj^ where a is the standard deviation associated with 

the noise process and the sum index 1=1, 2, • • • , m and the 

sum is over the values corresponding to the sample interval. 

This implies that 

Pf = 1 - P(x*/Ox) 

where P(z) is the standard normal cumulative distribution 
function. This relation is the basis for determining the 
threshold value x*. 

For the sample interval with signal energy, x is the value 

2 

of a normal random variable with a mean = 2 s^ where the sum 

index i = 1, 2, ■ ■ • , m and the sum is over the values 

corresponding to the sample interval. This implies that 

Pd = 1 - P(v* - d^) 

where v* = x*/Ox and d = S s^/a* = (l/m)-S/N. This relation 
is the basis for determining encounter detection probabilities 
for the first encounter model. The relation implies that for a 
specified false alarm probability pf the detection probability 
Pd is an increasing function of the signal to-noise ratio S/N. 
The Square Law Detector Statistic: The statistic for the second 

encounter model is a square law (energy) detector statistic that 
is defined by the sum 

X = 2 yi 

where the sum index i = 1, 2, • ■ • , m and the sum is over the 

values corresponding to the sample interval. For the second 
encounter model, only the characteristics of the noise are 
required to determine encounter detection probabilities. 
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is the 



For a sample interval without signal energy, x/a* 
value of a chi-square random variable with m degrees of 
freedom. This implies that 

Pf = 1 - P(x*/<^* U) 

where P(x*/a*|m) is the chi-square cumulative distribution 
function for a chi-square random variable with m degrees of 
freedom and where a is the standard deviation associated with 
the noise process. This relation is the basis for determining 
the threshold value x* . 

2 

For the sample interval with signal energy, x/a is the 

value of a noncentral chi-square random variable with m 

2 2 

degrees of freedom and noncentral parameter S s^/a where the 
sum index i = 1, 2, • • • , m and the sum is over the values 

corresponding to the sample interval. This implies that 

Pd = 1 ■ P(x*/cf*|ni. 2 Si/a*) 

JL 2 I 2 2, 

where P(x /a |m, E s^/a ) is the noncentral chi-square 
cumulative distribution function for a noncentral chi-square 

random variable with m degrees of freedom and noncentral 

2 2 , , . 
parameter S Sj^/a = (l/m)-S/N. This relation is the basis for 

determining encounter detection probabilities for the second 

encounter model. The relation implies that for a specified false 

alarm probability Pf, the detection probability p<^ is an 

increasing function of the signal-to-noise ratio S/N. This is 

made more evident by the following relation: 

P(x*/(j'|in, S s-/a') = 2{ (aVj!) -exp(-a) -PCX*/*?' I (m + 2- j) ] } 

where the sum index j = 0, 1, 2, ■ • • and the parameter 
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a 






(m/2) • (S/N) where the sum index 1=1, 2, 



= (1/2) -E Si/a = 

• • • , m. Note that the relation p^j > Pf is satisfied, 

P(x*/a* |m) < P[x*/a" | (m + 2 • j ) ] for j=0, 1, 2, ••• . 



since 
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Appendix 2. Program Probability Calculations 

The program evaluates the cumulative and inverse cumulative 
distribution functions using approximations described in 
Reference 2. These approximations are listed below. 

The Standard Normal Cumulative Distribution Function 
Approximation : 

P(z) = 1 - st-(bi + t- (b 2 + t- (b3 + t- (b4 + t-b5)))) 
where s = ( 1/ (2 • ir ) ^ ) • exp (-z V2 ) and t = 1/(1 + bQ-z). And 
where 

bo = .2316419, bi = .319381530, b2 = -.356563782, 
b4 = -1.821255978, b5 = 1.330274429, 
and z > 0. For z < 0, P(z) = 1 - P(|z|). 

The Inverse Standard Normal Cumulative Distribution Function 
Approximation : 

z(P) = t - (Cq + t-(Ci + tC 2 ))/(l + t(di + t- (d 2 + t-d 3 ))) 
where t = (ln(l/Q*))^ and Q = 1 - P. And where 
Cq = 2.515517, Cl = .802853, C 2 = .010328, 
di = 1.432788, d2 = .189269, d3 = .001308, 
and .5 < P < 1. For 0 < P < .5, z(P) = -z(l - P) . 

The Inverse Chi-Square Cumulative Distribution Function 
Approximation : 

v(p|m) = m-[l - 2/(9-m) + z • (2/ (9 m) ) I ] ^ 
where P(z) = P(v|m). In the program, the inverse standard 
normal cumulative distribution function approximation is used to 
determine z . 
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The Noncentral Chi-Square Cumulative Distribution Function 
Approximation : 

P(w|m,S s[/o") = P(z) 

where z= [2w/(l+b)]i -[2a/(l+b)-l]| with 

2 2 2 2 2 2 
a = m + S s^/a , b = (S s^/o )/(m + S s^/a ) and the sum index 

i=l, 2, •••,m. In the program, the standard normal 

cumulative distribution function approximation is used to 

determine P(z) . 
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Appendix 3. The Magnetic Signal 

The encounter models are defined by the following 
conditions: A submarine magnetic anomaly field is a magnetic 

dipole field that is superimposed on an earth magnetic field that 
is constant over an encounter region. A magnetometer magnetic 
signal value is the magnitude of the projection of a dipole 
magnetic field on the earth magnetic field at the location of the 
magnetometer. The basis for determining magnetic signal values 
in the program is an expression that can be developed as follows: 
In the right handed rectangular coordinate system that is shown 
in Figure 1, a magnetic dipole is at the origin, the xy-plane is 
the horizontal plane at a representative point in an encounter 
region, the positive y-axis is in the direction of magnetic 
north, the positive x-axis is in the direction of magnetic east 
and the positive z-axis is positive upward. In this rectangular 
coordinate system, the constant earth magnetic field can be 
expressed by = He(jcos § + k-sin $) where $ is a 

magnetic dip angle and Hg is a magnetic field magnitude that 
characterizes the earth field in an encounter region. 

In a spherical coordinate system with the origin at the 
magnetic dipole and the polar axis in the direction of the dipole 
moment, = (c-p/r^) • (2-r-cos 9 + 0-sin ©) is the magnetic 

field of the dipole at a point with spherical coordinates 
(r,r,9). In this expression p is the magnitude of the dipole 
moment and c is a constant whose value is determined by the 
choice of units. The magnetic signal for a magnetometer that is 
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z 




Figure 1. A unit vector in the direction of the earth field and a 
unit vector in the direction of the dipole field of a 
dipole at the origin are shown at a point that is a 
distance r from the dipole. The unit vector h is 
in the direction of the horizontal component of the 
dipole moment. 
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described by the encounter models is Hg = when the 

magnetometer is at the point (r,r,0) and Hg can be expressed in 
terms of the rectangular coordinate system as follows: First, 

let ro be a unit vector in the direction of the spherical 
coordinate system polar axis. Then the magnetic dipole field 
Hd = (c-p/r^) • (3-r-cos 0 - ro) 

since r = (ix + j-y + k-z)/r where r = (x* + y* + z*)^ , 

0 = [ (ro X r) x r]/sin 0 = [(r*ro)r - (r*r) • ro ]/sin 0, and the 
dot product r*ro = cos 0. The unit vector ro can be expressed 
in the rectangular coordinates by noting that ro = p/p and then 
expressing p in rectangular coordinates. To do this, let n 
be the depression angle of p from the xy-plane (the horizontal 
plane) with n positive downward and let a be the direction of 
p relative to magnetic north. Then the magnetic dipole moment 
p = p- (h- cos n - ksin n) in terms of the unit vector k and 
the unit vector h = i-sin a + j ■ cos n which has the direction 
of the horizontal component of p. With these results, the unit 
vector ro = i(cos n-sin a) + j • (cos n-cos a) - ksin n and 
“ (c-p/r^) • [ (3/r) -cos 0 (y cos $ - z-sin $) 

- (cos § • cos n-cos a + sin $-sin n) ] 
where cos © = (1/r) • (x-cos n-sin a + y-cos n-cos a - z-sin n) 
since cos © = r*rQ. As can be seen from this expression, for a 
constant dipole moment magnitude and direction and a constant 
earth field magnitude and direction, the magnetic signal is only 
a function of the rectangular coordinates of the location of the 
magnetometer relative to the dipole. In the encounter models. 



22 



both of these conditions are satisfied. However, 
n, a, Hg and $ to vary, the expression for Hg 
to more general encounter models. 



by allowing p, 
is applicable 
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Appendix 4 . The Anderson Formulation 

In the encounter models, the magnetic signal Hg at a 
sample point in a straight line encounter can be represented in a 
form described by Anderson in Reference 4. For convenience, the 
Anderson formulation is used in the program to determine values 
for Hg. It can be developed as follows: In Figure 2, the 

primed rectangular coordinate system is superimposed on the 
rectangular coordinate system of Figure 1 so that the origin is 
coincident with the origin of that system. A magnetometer is in 
a straight line encounter with a magnetic dipole that is located 
at the origin of the combined system. The combined system moves 
with the magnetic dipole with the x'-axis oriented so that it is 
parallel to and in the direction <f> of magnetometer's relative 
to the magnetic dipole and the z'-axis oriented so that it is 
directed toward and passes through the CPA on that track. Let 
1, m, and n be the direction cosines of the dipole moment p 
and 1^/ m^^ and n^ be the direction cosines of the earth 
magnetic field Hg in the primed system. Then the unit vector 
Tq = i'-l + j'n + k' -m. And, for points on the relative track, 
x' = s', y' =0 and z' = R where R is the slant range of the 
dipole at CPA and where s' is the algebraic distance of the 
magnetometer from CPA on the relative track. (It is negative for 
points before CPA and positive for points after CPA.) This 
implies that magnetometer position vector in the moving 
coordinate system is r = i'-(s'/r) + k' • (R/r) . 
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Magnetic 

East 



Figure 2 . 



The dipole is at the origin in the unprimed and primed 
coordinate systems. In the primed coordinate system, 
the CPA is at (0, 0, R) and, for a time t during 
the encounter, the magnetometer is at [s'(t), 0, R] . 
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From Appendix 3, 

Hd = (c-p/r^) • (3 • r- cos 0 - tq) • 

With r and tq expressed in terms of the primed unit vectors 
and cos 0 = r*rQ = l-(s'/i^) + h-(R/r), this becomes 

Hd = (cp/r3) • [ (3/r2) • (1-s' + n-R)(i'-s* + k'R) 

- (i' • 1 +j ' -m + k' -n) ] . 

Then, since Hg/Hg = i'-l^ + j'-m^ + k'-n^ and Hg = Hg-Hd/Hg, 
Hs = (c-p/r^) • [ (2 • 1 • li - m-mi - n-ni)-(s')* 

+ 3(nld + lni)s'-R + (2-n-ni - l-li - m-md)-R*]. 
The quantities 

A2 = 2-l-l]^ - m-m^ - n-n^ 

Ad ■= 3 • (n • Id + 1 • Hd) 

Aq = 2n-nd - I’ll “ m-md, 

are the Anderson coefficients. With r = [(s')* + R*]^ and 

P = s'/R, Hg can now be expressed as follows: 

HsO) = (C’P/r 3) -S Aj -Fj (P) 

where the Fj(p) = P^/(l + P*)^/^ for j = 0, 1, 2 are the 
Anderson functions. This is the Anderson formulation. 

To relate the Anderson formulation for Hg to the 
formulation for Hg in Appendix 3, first note that 

1 = ro*i' 
m = ro*j' 
n = ro*k' 

and 
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ll = (He/He) *i' 
mi = (He/He) -j' 
ni = (He/He) -k' . 

Then express Vq, (Hg/Hg) and the unit vectors i' , j' and k' 
in terms of the unit vectors i, j and k and take the 
indicated dot products. From Appendix 3, 

Tq = i-(cos n-sin a) + j • (cos n • cos a) - k-sin n 

and 

Hg/Hg = j • cos § + k-sin f>. 

To express i', j' and k' in terms of the unit vectors i, j 
and k, note that the unprimed coordinate system can be 
transformed to the primed coordinate system by two rotations that 
are defined as follows: First, rotate a coordinate system that 

is coincident with the unprimed coordinate system about the 
z-axis through the angle {<f> - ir/2) with positive angles clockwise 
(left hand rule) so that its x-axis is parallel to and in the 
direction of the relative track. Next, rotate this system about 
its x-axis through an angle 6 with positive angles 
counterclockwise (right hand rule) so that the positive z-axis 
passes through the CPA. The angle 6 is related to the 
vertical separation Zq of the magnetometer and the dipole and 
the algebraic encounter lateral range L that is positive if the 
dipole is to the left of the relative track. With these sign 
definitions: L = R-sin 6 and Zq = R-cos S. After the 

rotation, the auxiliary coordinate system is coincident with the 
primed coordinate system. 
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These transformations can be described in terms of matrix 
equations as follows; Let (x",y",z") be the coordinates of a 
point in the coordinate system that is coincident with the 
auxiliary coordinate system after the first rotation. Then the 
transformation from the unprimed coordinate system to this double 
primed coordinate system is described by the matrix equation 



X" 




sin <p cos 0 0 




X 


y" 


= 


-cos 0 sin 0 0 




y 


z" 




0 0 1 




z 



And the transformation from the double primed coordinate system 
to the primed coordinate system is described by the matrix 
equation 



X' 




10 0 




X" 


y 


= 


0 cos 5 sin 6 




y" 


z ' 




0 -sin 5 cos 6 




z" 



Taking the product of the rotation matrices in the indicated 
order yields the matrix equation 



X' 




sin 


0 




COS 


0 






0 






X 


y 


= 


-cos 


6 cos 


0 


cos 


6 


sin 


0 


sin 


5 




y 


z ' 




sin 


6 cos 


0 


-sin 


6 


sin 


0 


cos 


6 




z 



which defines the transformation from the unprimed to the primed 
coordinate system. The unprimed vector components of the unit 
vector i' can be found by transforming the coordinates (1,0,0) 
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in the primed system to their corresponding coordinates in the 
unprimed system with the inverse of this matrix and then 
repeating this process for (0,1,0) and (0,0,1) in order to find 
the unprimed unit vector components of j' and k' . However, 
since the inverse transformation matrix is the transpose of this 
matrix, the elements of the row that corresponds to a primed unit 
vector are the magnitudes of the unprimed vectors that are its 



components . 


Consequently: 


i' 


= i-sin (f> + j • cos <f> 


j' 


= -i-cos 5 cos <t> + j • cos S sin <p + k-sin 6 


k' 


= i-sin S cos <p - j-sin 6 sin <f> + k-cos 6. 


Then, taking 


the dot products between r, tq and these three 


unit vectors 


as indicate above gives; 


1 = 


= cos n-cos (0 - a) 


m = 


= cos 5 • cos n-sin {<p - a) - sin 6 -sin n 


n = 

and 


= -sin 6 • cos n-sin (<p - a) - cos 6-sin n 


ll 


= cos § - cos 0 




= cos 6 - cos ^>-sin 0 - sin 6 -sin $ 


ni 


= -sin 6 - cos $-sin 0 - cos 6-sin ^> . 



These are the relations that are used in the program to determine 
values for the Anderson coefficients. 
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Appendix 5. The Encounter Equations of Motion 

In the double primed coordinate system that is defined in 
Appendix 4, the equations of motion of a magnetometer relative to 
a submarine (dipole) are: 
x”(t) = s' (t) 
y"(t) = -L 
z"(t) = Zq. 

where L is the algebraic encounter lateral range that is 
defined in Appendix 4, zq is the vertical separation between the 
magnetometer and the submarine and s'(t) is the distance of the 
magnetometer from the CPA on the relative track. With w the 
speed of the magnetometer relative to the submarine and t a 
relative time parameter, s'(t) = w-t . These equations can be 
considered to be the ones used in the program to describe the 
motion of a magnetometer relative to a submarine. There, t is 
determined by t = [j - (m-l)/2]St where the index j = 1, 2, 

• • • , m and 5t, a time step, is the time between samples. Note, 
when t = 0, the magnetometer is at the CPA. 

In the coordinate system of Figure 1 in a straight line 
encounter as defined in the encounter models, the equations of 
motion of a magnetometer relative to a submarine can be written 
as follows: 

x(t) = s' (t) -sin 0 + L-cos 0 
y(t) = s' (t) cos 0 - Lsin 0 

Z(t) = Zq, 
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since the transformation from the double primed coordinate system 
to the primed coordinate system is determined by the matrix 
equation 



X 




sin 0 -cos 0 0 




X" 


y 




cos 0 sin 0 0 




y" 


z 




0 0 1 




z” 



X" = s'(t) and y” = -L. The above equations and the expression 
for Hg in Appendix 3 could have been used in the program to 
evaluate the magnetic signal. In particular, with these 
equations of motion and with the two relations L = Rsin 5 and 
Zq = R-cos 5, the expression for Hg in Appendix 3 can be 
written in terms of f2, a, §, <p, S and R so that it is 
identical in appearance to the Anderson formulation for Hg in 
terms of these quantities. The definition of 6 in Appendix 4 
in terms of a counterclockwise rotation results in a definition 
of the algebraic lateral range that is consistent with some that 
have been used elsewhere. 

In the program, the relation w = v - u is the basis for 
determining the relative speed w. In this relation, v is the 
velocity of the magnetometer, u is the velocity of the submarine 
(dipole) and w is the velocity of the magnetometer relative to 
the submarine. This relation implies the following equations: 

Wx = Vx - Uj^, Wy = Vy - Uy and w^ = v^ - U 2 where the 
coordinates x, y and z refer to a fixed coordinate system 
with the same orientation as that of Figure 1. In the encounters 
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of the models, Vy = vsin a, Vy = v • cos a and V 2 = 0 where a 
is the magnetic course and v is speed of the magnetometer. 

And, Ujj = u-sin p, Uy = u cos p and U 2 = 0 where p is the 
magnetic course and u is the speed of the submarine. The 
relative magnetic course <p and the relative speed of the 
magnetometer are defined by = w-sin 0 and Wy = w-cos <p. In 

the program, 0 and w are determined with a rectangular to 
polar conversion routine where w = ( w^ + Wy) ^ and where <p is 
determined by sin”^(Wjj/w) and cos"^(Wy/w). 
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Appendix 6. The Sxibmarine Magnetic Dipole 

In the encounter models, the magnitude and direction of a 
submarine's dipole moment p are determined by first determining 
its components in the rectangular coordinate system of Figure 1. 
In that coordinate system, 

Px = (PLP + PLI) *sin P + (Ptp + Pti) ' cos p 
Py = (Plp Pli) 'COS p - (Ptp + Pti) ’Sin p 
Pz ^ ~(PVP Pvi) 

where p is the submarine's magnetic course, PLp, Ptp pyp 

are the permanent and Plj; Pti Pvi induced 

longitudinal, transverse and vertical magnetic dipole moments of 
the submarine. These relations are based on the following sign 
convention: 



PL 


is 


positive 


when 


PL 


is 


directed 


from 


stern to bow. 


pt 


is 


positive 


when 


Pt 


is 


directed 


from 


port to starboard. 


Pv 


is 


positive 


when 


Pv 


is 


directed 


downward. 



The permanent dipole moments are input parameters in the program 
and the induced dipole moments are determined in a way that is 
similar to one that is described in Reference 1. In the 
encounter models, a submarine is a ferromagnetic prolate 
ellipsoid with the major axis the submarine's longitudinal axis 
and the equal minor axes the submarine's transverse and vertical 
axes. And, the induced dipole moments are: 

Pli = ^L'J^eL 

Pti = kT'HgT 

Pvi = ky-Hev 
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where Hgip and are the vector components of the earth 

magnetic field in the rectangular coordinate system defined by 
the ellipsoid axes and the submarine magnetic moment sign 
convention. The magnitude of these vector components are: 

HgL = Hg-cos $ cos p 
Hgip = -Hg-cos $ sin p 
HgY “ Hg • sin 

The earth magnetic field dip angle $ and the earth magnetic 
field magnitude Hg can each be chosen to be input parameters in 
the program or the can be computed by the program as described in 
Appendix 7. (Values for $ and Hg can be found from magnetic 
data charts, for example, see References 5 and References 6.) 

By using the above relations, 

Px = Hg-cos (Kl - kij.) - sin p cos p + (p^p sin p + p.pp-cos P) 

Py = Hg-cos § (1cl-cos*P + k.p sin*P) + (PLp-cos p - pipp-sin P) 
pz = -(Hg-kv-sin § + pvp) 

In the encounter models, the values of the permeability 
coefficients kL, krp and ky are related to submarine 
displacement by the following relations: 
kL = fL-W 
kiji — fip • w 

ky “ f Y • W 

where W is the submarine displacement in tons and fj_^, fq, and 
fY are permeability factors that are determined by a submarine's 
magnetic characteristics. In the program, the units of Hg are 
gamma, the units of p are in oersted-centimeter^, the units of 
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k are oersted-centimeter^/gamma and the units of f are 
oersted-centimeter^/gamma-ton . If the units of Hg were 
gamma, but the units of p were gamma-foot^, then the units of 
k would be foot^ and the units of f would be foot^/ton. To 
convert p in gamma-foot^ to oersted-centimeter^ or k in 
foot^ to oersted-centimeter^/gamma or f in foot^/ton to 
oersted-centimeter^/gamma-ton, divide by 3.53. (The program 
default values for fj^, fip and fy are values from Reference 1 
in foot^/ton that have been divided by 3.53 to give values in 
oersted-centimeter^/gamma-ton. ) To convert p in weber-meter to 
oersted-centimeter^, multiply by [ 1/ (4 ■ ti) ] • 10^*^ and to convert 
p in ampere-centimeter^ to oersted-centimeter^, multiply by 
1 - 10 ^ . 
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Appendix 7. The Earth Magnetic Field and Dip Angle 

An auxiliary magnetic field model is described in this 
appendix. The model is the basis for a default choice for either 
the value of the earth magnetic field magnitude parameter Hg or 
the earth magnetic field dip angle parameter Relative to 

encounter model accuracy, the default values should be adequate 
in most cases. 

In the model, the earth magnetic field is generated by a 
magnetic dipole that is located at the earth's center and the 
earth is a nonmagnetic sphere of radius r@. With pg the 
magnitude of the dipole moment, the magnitude of the earth field 
at a point is 

He = (He-He)^ = (c • Pg/r^ ) • ( 3 • cos *6 + 1 )^ 
where 6 is the polar angle of the point in a spherical 
coordinate system and c is a constant whose value is determined 
by the choice of units. The dipole moment is coincident with and 
in the direction of the polar axis which is directed toward the 
earth's southern hemisphere. In this coordinate system, at any 
point on the surface of the earth: 

He = HeQ-(3-cos*0 + 1)^ 

where H^q is the value of Hg at the magnetic equator which is 
defined by the points on the earth where 0 = 90°. In terms of 
the dip angle, at any point on the earth's surface: 

Hg = 2 -HgQ- (3 -cos*# + 1)“^. 

This expression can be obtained by noting that # can be 
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defined in terms of the or the r and 0 components of Hg 
the spherical coordinates as follows: 



in 



sin § = -Hgj-/Hg = -2-Hgo-(cos e)/Hg 

and 

cos § = He©/Hg = Hgo- (sin 6)/He. 

Based on these relations, 

[(sin $)/2]’ + cos"$ = (Hgo) * 

and 

Hg = HgQ ■ [ (sin ' ) /4 + cos*$]“2 = 2 HgQ- (3 • cos*§ + 1)“^ 
The dip angle is determined from a magnetic latitude for the 
encounter region. 

The magnetic latitude and longitude of the representative 
point of an encounter region can be defined in terms of its 
geographic latitude and longitude by using the following 
transformations: First, convert the latitude and longitude of 

the point to rectangular coordinates in a right-handed coordinate 
system whose origin is at the center of a spherical earth, whose 
positive z-axis passes through its north geographic pole and 
whose positive y-axis passes through the point on its equator 
with latitude 0° and longitude 0°. Next, rotate a coordinate 
system that is coincident with this system in a clockwise 
direction (left hand rule) about its z-axis so that its positive 
y-axis passes through the point with latitude 0° and with 
longitude equal to that of the north magnetic pole. Then, rotate 
the system in a clockwise direction about its x-axis so that its 
positive z-axis passes through the north magnetic pole. Next, 



37 



transform the rectangular coordinates of the representative point 
in this system to its coordinates in the spherical coordinate 
system that is associated with it. Then, with n the polar 
angle of the representative point in this system, the magnetic 
latitude of the representative point is Lm “ Since 

tan = -2 • cot © from above and 9 = 180° ~ fl, the dip angle $ 
is given by the following relation: 

$ = tan“^(2-tan L^) . 

In the program, the transformations described above are 
accomplished in part by a rectangular to polar conversion 
routine. In particular, by rotating the final rectangular 
coordinate system about its z-axis so that the x-coordinate is 
zero, the polar angle of the polar coordinates of the 
representative point in the resulting yz-plane determines Ljjj. 

Some values of $ and Hg that are listed in Table 1 were 
generated using the program. The latitude and longitude of the 
magnetic pole that are in the program are 76° N and 100° W. 
Values of $ from Reference 6 and of Hg from Reference 7 are 
also listed in Table 1. A comparison of the program values with 
these values gives an indication of the errors inherent in the 
procedure. 
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Encounter 

Latitude 

60°N 

60°N 

60°N 

60°N 

30°N 

30“N 

30°N 

30°N 

30°N 

30°N 

30“N 

0 ° 

0 ° 

0 ° 

0 ° 

30°S 

30°S 

30°S 

60°S 



Encounter 

Longitude 

180°W 

150°W 

90°W 

30°W 

150°W 

90”Vl 

60°W 

30°E 

60°E 

90°E 

150°E 

90°W 

30°W 

0 “ 

60°E 

180°W 

90°W 

90°E 

30°W 



Program 
Dip $ 

74 ° 
78° 
82° 
75° 
58° 
62° 
59° 
37° 
31° 
30° 

42 ° 
26° 

9° 

-48° 

-25° 

-45° 

-30° 

-62° 

-69° 



Chart #30 
Dip 4 

70° 

73 ° 

84° 

74 ° 

50° 

60° 

59° 

43° 

45° 

44° 

40° 

18° 

-5° 

-25° 

- 20 ° 

-55° 

-34° 

- 66 ° 

-56° 



Program 
Field He 

. 63 
. 66 
. 68 
. 64 
. 51 
. 55 
. 53 
.41 
.39 
. 39 
.43 
.38 
.35 
.35 
. 38 
.44 
.39 
.55 
. 60 



Chart #39 
Field He 

.53 
. 55 
.61 
.52 
.41 
. 51 
.48 
. 42 
.46 
. 44 
.40 
. 32 
. 28 
.31 
. 37 
. 50 
.31 
. 52 
.34 



Table 1. Some corresponding values for the encounter dip angle 
§ and the earth magnetic field intensity He in 
gamma. Chart #33 refers to Reference 5 and Chart #39 
refers to Reference 6. 



A small difference between a chart value and a computed value for 
the dip angle $ at some location does not imply a small 
difference between the chart value and the computed value for the 
earth magnetic field intensity at that location. For example, 
note the values for 60°N latitude and 30°W longitude. 
(Positive values of $ indicate the inclination or dip below the 
horizontal of the north seeking end of a dip needle. Negative 
values indicate the inclination below the horizontal of the south 
seeking end.) 
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An extension of the above procedure can be made for finding 
the magnetic variation at a location. However, the values 
generated by using the procedure are generally unsatisfactory. 
Magnetic variation values are charted in Reference 7. 
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Appendix 8 . An Alternative Encounter Model 

The alternative model that is referred to in Section III is 
described in more detail in this appendix. In the model, the 
detection range is the slant range R at the CPA for an 
encounter with a specified detection probability (usually .5) is 
defined by: 

R = [c-p/Hs]V3 

where c and p are defined in Appendix 3 and Hg represents 
a minimum detectable average magnetic signal that is defined by: 
Hg = (ORF) -Nm 

where ORF is a signal-to-noise ratio called the operator 
recognition factor and represents the magnetic noise. 

Combining these two relations gives: 

R = {c-p/[ (ORF) -N^] } V3. 

The value for ORF depends both on the specified detection 
probability and on a specified or implied false alarm 
probability. 

The Anderson formulation is consistent with these relations 
in an approximate sense if the average magnetic signal H is 
defined as a root mean square value such that H = (cp/R^) -k 
where k is an encounter parameter defined by: 
k = {Z [Z Aj -Fj (Pi) 3 

with the first sum index i = 1, 2, • • • , m and the second sum 

index j = 0, 1, 2. For a particular encounter geometry, k is 
constant and this suggests that the two encounter models could be 
used to determine an average value for k for an encounter 
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region based on average submarine magnetic characteristics. 

Values for both k and H are generated by the program and such 
values can give an indication of the magnitude of the differences 
in detection range estimates that are based on this model and 
those that are based on either of the other two encounter models. 
A more detailed comparison of these encounter models is described 
in Reference 8. 
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Appendix 9. An Example of the Program Output 

The program is designed to generate the following 
cpaantities : encounter parameter values, lateral range function 
values for the crosscorrelation encounter model and for the 
square law encounter model, average magnetic signal values, slant 
range at CPA values, encounter parameter values and magnetic 
signal values. These values can be saved as a program data file 
and/or they can be printed. 

An example of the program's printed output is listed in 
Table 2, Table 3 and Table 4. Figure 3 is a plot of the lateral 
range function values (a lateral range curve) for the square law 
detector that are listed in Table 3. Figure 4 is a plot of the 
magnetic signal values that are listed in Table 4. 

A second example of the program's printed output is listed 
in Table 5, Table 6 and Table 7. Figure 5 is a plot of the 
lateral range function values for the square law detector that 
are listed in Table 6. Figure 6 is a plot of the magnetic signal 
values that are listed in Table 7. 
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progran file name 
program data file name 
magnetic data file name 
processing data file name 
kinematic data file name 



MAD. BAS 
data .mad 
data .mag 
data. pro 
data. kin 



combined magnetic, processing & kinematic data file name data.ropk 



encounter latitude (decimal degrees) 
encounter longitude (decimal degrees) 
encounter variation (decimal degrees) 
encounter dip angle (decimal degrees) 
encounter magnetic field intensity (gamma) 
permanent longitudinal moment (oersted-cm3 ) 
permanent transverse moment (oersted-cm3 ) 
permanent vertical moment (oersted-cm3 ) 
target displacement (tons) 

target longitudinal permeability coefficient 
target transverse permeability coefficient 
target vertical permeability coefficient 
target longitudinal permeability factor 
target transverse permeability factor 
target vertical permeability factor 
sampling period (seconds) 
integration time (seconds) 
adjusted integration time (seconds) 
number of samples per encounter 
magnetometer course (decimal degrees) 
magnetometer speed (knots) 
magnetometer altitude (meters) 
target course (decimal degrees) 
target speed (knots) 
target depth (meters) 

magnetometer relative course (decimal degrees) 

magnetometer relative speed (knots) 

magnetometer-target vertical separation (meters) 

target induced longitudinal dipole moment (oersted-cm3 ) 

target induced transverse dipole moment (oersted-cm3 ) 

target induced vertical dipole moment (oersted-cm3 ) 

magnetic dipole moment (oersted-cm3 ) 

dipole moment azimuth (decimal degrees) 

dipole moment depression angle (decimal degrees) 

distance between samples on the relative track (meters) 

false alarm rate (false alarms per hour) 

false alarm probability 

magnetic noise (gamma) 

maximum lateral range (meters) 

lateral range step (meters) 

number of lateral range function values 



45 

-45 

-25 

68.41608 

59035.05 
0 

0 

0 

1800 

13140 

2880 

2880 

7.3 

1.6 

1.6 

.5 

20 

20.5 
41 
290 
180 
100 
20 
10 
100 

286.8202 

180.2776 

200 

2.017795E+08 

-4.422564E+07 

1.58099E-I-08 

2.601272E+08 

7.637512 

37.42885 

46.3714 

3 

1. 666667E-02 

.35 

500 

20 

51 



Table 2. 



An example of an encounter parameter values printout. 
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data.r.ad lateral range function values 



L 


p(cc) 


p(sl) 


H 


R 


k 


meters 






gamma 


meters 




-500 


.2155959 


. 0255824 


. 184247 


538.5165 


2.818488 


-480 


.2721803 


. 029692 


.21034 


520 


2.879916 


-460 


. 3464673 


3.582633E-02 


. 2410144 


501.5974 


2.943744 


-440 


.4415205 


4 . 529992E-02 


.2785504 


483 . 3218 


3 . 009938 


-420 


. 5573696 


6 . 047795E-02 


. 3231911 


465. 1881 


3 . 078418 


-400 


. 6872897 


. 0856887 


. 3764825 


447.2136 


3 . 149025 


-380 


.8146549 


. 1287804 


.4404623 


429.4182 


3.221497 


-360 


.9159235 


.2030898 


.5177304 


411.8252 


3 . 29544 


-340 


.9745625 


. 3272475 


.6108596 


394.4617 


3 . 370274 


-320 


.9958122 


.5141411 


.7233319 


377. 3593 


3 .445185 


-300 


. 9997211 


.7388175 


.8597726 


360.5551 


3.519032 


-280 


.9999951 


.9177597 


1.025158 


344.093 


3 . 590261 


-260 


1 


. 9896971 


1.224529 


328 . 0244 


3.656775 


-240 


1 


.999714 


1.463644 


312.41 


3 .71578 


-220 


1 


.9999992 


1.74981 


297. 3214 


3 .763606 


-200 


1 


1 


2.085394 


282.8427 


3.795517 


-180 


1 


1 


2.470496 


269.0725 


3.805514 


-160 


1 


1 


2.901352 


256.125 


3.786234 


-140 


1 


1 


3.359555 


244 . 1311 


3.72898 


-120 


1 


1 


3 .819509 


233.2381 


3.624058 


-100 


1 


1 


4 . 239317 


223 . 6068 


3.461596 


-80 


1 


1 


4.532758 


215.4066 


3 .232958 


-60 


1 


1 


4.618943 


208.8061 


2.932828 


-40 


1 


1 


4.44342 


203.9608 


2.5617 


-20 


1 


1 


3.97739 


200.9975 


2.128285 


0 


1 


1 


3 . 332139 


200 


1.651697 


20 


1 


.9999993 


2 . 540665 


200.9975 


1.165581 


40 


.9999931 


.908713 


1.76834 


203.9608 


.7394171 


60 


.9933111 


.467368 


1.110188 


208 . 8061 


. 5637183 


80 


. 9998733 


. 7874625 


.9039734 


215.4066 


.7782836 


100 


.9999999 


.9856169 


1.086971 


223 . 6068 


1. 134662 


120 


1 


.9992399 


1.254863 


233.2381 


1 . 490748 


140 


1 


. 9998478 


1 . 308583 


244 . 1311 


1 .811734 


160 


1 


. 9998431 


1 . 29482 


256.125 


2 . 089946 


180 


1 


.9994995 


1 . 230946 


269.0725 


2.326072 


200 


1 


. 9974561 


1.131703 


282.8427 


2.523877 


220 


1 


.9872744 


1 . 017684 


297.3214 


2 . 688154 


240 


.999999 


. 9510013 


. 9018685 


312.41 


2 . 823785 


260 


.9999747 


.8646194 


.7917129 


328.0244 


2.935307 


280 


.9996544 


.7242736 


.6909791 


344 . 093 


3.026728 


300 


.9973592 


. 5582789 


.6011183 


360.5551 


3.10149 


320 


. 9876337 


.4037092 


.5222269 


377.3593 


3.162495 


340 


.9608865 


. 2816961 


.4536729 


394 .4617 


3.212155 


360 


.9086407 


. 1948741 


. 3944811 


411.8252 


3.252463 


380 


. 8301703 


. 1365131 


.3443362 


429.4182 


3.285058 


400 


. 7333605 


9.820431E-02 


. 3014853 


447.2136 


3.311281 


420 


. 6299387 


. 0731037 


. 2645045 


465.1881 


3.332228 


440 


. 5302704 


5 . 646636E-02 


.2325951 


483.3218 


3.348799 


460 


. 4408869 


4 . 522933E-02 


.205043 


501.5974 


3.361725 


480 


. 3645446 


3.747106E-02 


. 1812233 


520 


3 . 371603 


500 


. 3013791 


3. 199235E-02 


. 1605954 


538.5165 


3 . 378925 



Table 3. An example of a lateral range function values printout. 

The heading for the crosscorrelation values is p(cc) 
and the heading for the square law values is p(sl) • 
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Figure 3. A plot of the square law lateral range function values 
that are listed in Table 3. The horizontal axis is 
encounter horizontal range at CPA in meters. The 
vertical axis is encounter detection probability. 
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data.n^ad riagnetic signal values 
for a lateral range of 40 meters 



relative CPA distance 


magnetic signal 


s in meters 


Hs in gamma 


-927.4282 


-1.619177E-02 


-881. 0568 


-1.834315E-02 


-834 . 6854 


-2.086144E-02 


-788.314 


-2 . 381637E-02 


-741.9426 


-2 .728804E-02 


-695.5712 


-.0313641 


-649 . 1998 


-3 . 613118E-02 


-602.8283 


-4.165412E-02 


-556.4569 


-4 .792992E-02 


-510.0855 


-5. 478963E-02 


-463.7141 


-6. 169327E-02 


-417.3427 


-6. 730653E-02 


-370.9713 


-.0686387 


-324 . 5999 


-5.933568E-02 


-278 . 2285 


-2 . 652314E-02 


-231.8571 


. 0540608 


-185.4856 


.2215785 


-139. 1142 


. 5163918 


-92.74282 


.9104893 


-46.37141 


1.184586 


0 


1.003207 


46.37141 


. 3837067 


92 . 74282 


-.2252113 


139.1142 


-.5323231 


185.4856 


-.5837545 


231.8571 


-.5168925 


278 . 2285 


-.4202569 


324 . 5999 


-.3304064 


370.9713 


-.2572838 


417.3427 


-. 2007385 


463 .7141 


-. 1577988 


510.0855 


-.125285 


556.4569 


-. 1005532 


602.8283 


-8 . 158605E-02 


649.1998 


-6. 689598E-02 


695.5712 


-5. 539922E-02 


741. 9426 


-4 . 630719E-02 


788.314 


-3.904362E-02 


834 . 6854 


-3.318428E-02 


881.0568 


-.0284142 


927.4282 


-2 . 449729E-02 



Table 4. An example of a magnetic signal values printout. 
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Figure 4. A plot of the magnetic signal values listed in Table 4 
The horizontal axis is relative CPA distance in meters 
The vertical axis is magnetic signal in gamma. 
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program file name 
program data file name 
magnetic data file name 
processing data file name 
kinematic data file name 

combined magnetic, processing & kinematic data file name 

encounter latitude (decimal degrees) 

encounter longitude (decimal degrees) 

encounter variation (decimal degrees) 

encounter dip angle (decimal degrees) 

encounter magnetic field intensity (gamma) 

permanent longitudinal moment (oersted-cm3 ) 

permanent transverse moment (oersted-cm3 ) 

permanent vertical moment (oersted-cm3 ) 

target displacement (tons) 

target longitudinal permeability coefficient 
target transverse permeability coefficient 
target vertical permeability coefficient 
target longitudinal permeability factor 
target transverse permeability factor 
target vertical permeability factor 
sampling period (seconds) 
integration time (seconds) 
adjusted integration time (seconds) 
number of samples per encounter 
magnetometer course (decimal degrees) 
magnetometer speed (knots) 
magnetometer altitude (meters) 
target course (decimal degrees) 
target speed ()cnots) 
target depth (meters) 

magnetometer relative course (decimal degrees) 

magnetometer relative speed (knots) 

magnetometer-target vertical separation (meters) 

target induced longitudinal dipole moment (oersted-cm3 ) 

target induced transverse dipole moment (oersted-cm3 ) 

target induced vertical dipole moment (oersted-cm3 ) 

magnetic dipole moment (oersted-cm3 ) 

dipole moment azimuth (decimal degrees) 

dipole moment depression angle (decimal degrees) 

distance between samples on the relative track (meters) 

false alarm rate (false alarms per hour) 

false alarm probability 

magnetic noise (gamma) 

maximum lateral range (meters) 

lateral range step (meters) 

number of lateral range function values 



MAD. BAS 
datal , mad 
data . mag 
data .pro 
datal . kin 
datal . mpk 
45 
-45 
-25 

68.41608 

59035.05 
0 

0 

0 

1800 

13140 

2880 

2880 

7.3 

1.6 

1.6 

.5 

20 

20.5 
41 
290 
180 
100 
273 
10 
100 

290.9828 

170.462 

200 

1. 339682E+08 

5.522353E-^07 

1.58099E+08 

2 . 144584E+08 

295.4022 

47.49356 

43.84662 

3 

1. 666667E-02 
.35 
500 
20 
51 



Table 5. An example of an encounter parameter values printout. 
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datal.jnad lateral range function values 



L 


p(cc) 


p(sl) 


H 


R 


k 


meters 




gamma 


meters 




-500 


9.718663E-02 


1. 815607E-02 


.1162113 


538.5165 


2.117317 


-480 


. 1097076 


1.887764E-02 


.1307452 


520 


2.066198 


-460 


.1252767 


1. 979502E-02 


. 1476683 


501.5974 


2.017377 


-440 


. 1449347 


2 . 098504E-02 


.1677343 


483 . 3218 


1.971743 


-420 


. 1701671 


2.256657E-02 


.1915806 


465.1881 


1.930378 


-400 


.2031047 


2 . 473001E-02 


.2196898 


447.2136 


1.89458 


-380 


.2467793 


2 .779396E-02 


.2543952 


429.4182 


1.865873 


-360 


.3053632 


3 . 231674E-02 


.2951223 


411.8252 


1.845999 


-340 


. 3841007 


3 .932704E-02 


.3461031 


394.4617 


1.836884 


-320 


.4881344 


5.082678E-02 


.4085031 


377.3593 


1.840555 


-300 


.6185669 


7. 091362E-02 


.4830141 


360.5551 


1.859016 


-280 


.7643011 


.1082717 


.5793262 


344 . 093 


1.89408 


-260 


.8946846 


. 1810074 


. 6993423 


328.0244 


1.947158 


-240 


.9731281 


. 3214172 


.8450639 


312.41 


2.019043 


-220 


.9973629 


.5584107 


1.022772 


297.3214 


2.109683 


-200 


.9999472 


.8329318 


1.262045 


282.8427 


2.217957 


-180 


.9999999 


.9794214 


1.555326 


269.0725 


2.341425 


-160 


1 


.9996887 


1.906387 


256.125 


2.476081 


-140 


1 


. 9999999 


2.315586 


244 . 1311 


2.616057 


-120 


1 


1 


2.774995 


233.2381 


2.753399 


-100 


1 


1 


3.289594 


223.6068 


2.878014 


-80 


1 


1 


3.848771 


215.4066 


2.978009 


-60 


1 


1 


4.356607 


208.8061 


3.040665 


-40 


1 


1 


4.738054 


203 . 9608 


3.054164 


-20 


1 


1 


4.921432 


200.9975 


3.009893 


0 


1 


1 


4.862248 


200 


2.904766 


20 


1 


1 


4.561756 


200.9975 


2 . 742697 


40 


1 


1 


4.06876 


203.9608 


2.534609 


60 


1 


1 


3.462592 


208.8061 


2.296925 


80 


1 


1 


2.876776 


215.4066 


2.049263 


100 


1 


. 9999861 


2.341757 


223.6068 


1.812213 


120 


1 


.9908414 


1.862014 


233.2381 


1.60562 


140 


.9999632 


.8492265 


1.463262 


244 . 1311 


1.446881 


160 


.9973499 


.55795 


1.150695 


256.125 


1.348197 


180 


.9773154 


.3394268 


.9187382 


269.0725 


1.31296 


200 


.9311967 


. 2232209 


.7327372 


282.8427 


1.334314 


220 


.8711751 


. 1618173 


.5979505 


297.3214 


1.39856 


240 


. 8071772 


. 1253278 


.4948332 


312.41 


1.49089 


260 


.7414287 


. 1006817 


.4163041 


328.0244 


1.599193 


280 


.6739562 


8.250537E-02 


. 3566572 


344.093 


1.714996 


300 


. 6054811 


.0684997 


.3130742 


360.5551 


1.832886 


320 


.5376865 


5.753178E-02 


.2806611 


377.3593 


1.949611 


340 


.4726234 


4.891159E-02 


.2530364 


394.4617 


2.063321 


360 


.4121077 


4.213673E-02 


.2280127 


411.8252 


2.173046 


380 


.3573765 


3.681085E-02 


.2054994 


429.4182 


2.27835 


400 


. 3090144 


3. 261634E-02 


.1853368 


447.2136 


2.379115 


420 


.2670518 


.0293013 


. 1673323 


465.1881 


2.475412 


440 


.2311289 


2.666843E-02 


.151282 


483.3218 


2.567413 


460 


.2006615 


2.456539E-02 


.1369845 


501.5974 


2.655341 


480 


. 1749714 


. 0228749 


. 12437 


520 


2.739442 


500 


. 1533755 


2 . 150718E-02 


.1132715 


538.5165 


2.819965 



Table 6. 



An example of a 
The heading for 
and the heading 



lateral range function values printout, 
the crosscorrelation values is p(cc) 
for the square law values is p(sl). 
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Figure 5. A plot of the square law lateral range function values 
that are listed in Table 6. The horizontal axis is 
encounter horizontal range at CPA in meters. The 
vertical axis is encounter detection probability. 
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datal.mad magnetic signal values 
for a lateral range of 40 meters 

relative CPA distance magnetic signal 

s in meters Hs in gamma 



-876.9324 

-833.0858 

-789.2391 

-745.3925 

-701.5459 

-657.6993 

-613.8527 

-570.006 

-526.1594 

-482.3128 

-438.4662 

-394.6196 

-350.7729 

-306.9263 

-263.0797 

-219.2331 

-175.3865 

-131.5399 

-87.69324 

-43.84662 

0 

43.84662 

87.69324 

131.5399 

175.3865 

219.2331 

263.0797 

306.9263 

350.7729 

394.6196 

438.4662 

482.3128 

526.1594 

570.006 

613.8527 

657.6993 

701.5459 

745.3925 

789.2391 

833.0858 

876.9324 



6.869243E-03 
8.988515E-03 
1. 181241E-02 
1.561336E-02 
2.078687E-02 
2 .791575E-02 
3.787262E-02 
5.198577E-02 
7.231193E-02 
.1020894 
. 1464967 
.2139163 
.317984 
.480659 
.735856 
1.130141 
1.707647 
2.448238 
3.132554 
3.263226 
2.439828 
1.023829 
-.1487533 
-.7045798 
-.8055338 
-.7078842 
-.5631492 
-.4307948 
-.3256765 
-.2465582 
-.1881302 
-.1451051 
-.1132614 
-8.948089E-02 
-7. 152954E-02 
-5.782279E-02 
-4.723642E-02 
-3.896839E-02 
-3.244173E-02 
-2.723733E-02 
-2.304761E-02 



Table 7. 



An example of a magnetic 



signal 



values printout. 
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3,5 




Figure 6. A plot of the magnetic signal values listed in Table 
7. The horizontal axis is relative CPA distance in 
meters. The vertical axis is magnetic signal in 
gamma . 
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Appendix 10. A Program Listing 
10 CIS : CLEAR 

20 EFJINT "Magnetic Anomaly Detection (MAD) Lateral Range Function Program" 

30 DIMX0(70), R0(70), H(70) , HS(70, 150), K(70) , HXX:(70) , PDSL(70) 

40 PI = 3.141592654#: CON = 2 * PI / 360: KDN = 1852 / 3600: N$ = "MAD.EAS" 

50 QO = .2316419: Q1 = .31938153#: Q2 = -.356563782#: Q3 = 1.781477937# 

60 Q4 = -1.821255978#: Q5 = 1.330274429# 

70 II = 2.515517: 12 = .802853: 13 = .010328: 14 = 1.432788: 15 = .189269: 16 
= .001308 
80 PRINT 

90 INRJT "generate data or print a program data file (g/p) " ; A$ 

100 IF A$ = "G" OR A$ = "g" THEN GOTO 120 

110 IF A$ = "P" OR A$ = "p" THEN GOTO 3070 EISE GOTO 90 

120 PRINT : A$ = "a" 

130 INRJT "magnetic, processing & kinematic data entry by combined file 
(y/n)"; A$ 

140 IF A$ = "Y" OR A$ = "y" IHEN GOTO 2120 

150 IF A$ = "N" OR A$ = "n" THEN GOTO 160 ELSE GOTO 130 

160 INRJT "magnetic data entry by file/keitoard (f/k)"; A$ 

170 IF A$ = "F" OR A$ = "f" THEN GOTO 780 

180 IF A$ = "K" OR A$ = "k" IHEN GOTO 190 ELSE GOTO 160 

190 INRJT "latitude in decimal degrees (N +)"; LAT: lATR = lAT * OCa4 

200 INRJT "longitude in decimal degrees (E +)"; ING: INGR = -ING * OON 

210 INRJT "variation in decimal degrees (E +)"; DEC: DECR = DEC * CON 

220 A$ = "a" 

230 INRJT "input dip angle (y/n)"; A$ 

240 IF A$ = "Y" OR A$ = "y" IHEN GOTO 260 

250 IF A$ = "N" OR A$ = "n" IHEN GOTO 290 ELSE GOTO 230 

260 INRJT "dip angle in decimal degrees (north magnetic hemisphere +)"; DIP 
270 DIFR = DIP * OON 
280 GOTO 390 

290 lATPR = 76 * CON: INGFR = 100 * OON 

300 D = SIN(INGR - INGPR) * OOS(IATR) : E = COS (INGR - INGPR) * COS(IATR) : F = 
SIN (LATR) 

310 X = E: Y = F 
320 GOSUB 4170 

330 T = T - (90 * cm - lATIR) : E = R * SIN(T) : F = R * OOS(T) 

340 X = E: Y = D 
350 GOSUB 4170 
360 X = F: Y = R 
370 GOSUB 4170 

380 DIPR = AIN(2 * TAN(T)) : DIP = DIFR / Cm 
390 A$ = "a" 

400 INFUT "input encounter magnetic field intensity (y/n)"; A$ 

410 IF A$ = "Y" OR A$ = "y" THEN GOTO 430 
420 IF A$ = "N" OR A$ = "n" IHEN GOTO 450 ELSE GOTO 400 
430 INFUT "encounter magnetic field intensity in gamma"; HE 
440 GOTO 460 

450 HE = 70000! / SQR(3 * COS (DIFR) * COS(DIFR) + 1) 

460 PLM = 0: PIM = 0: PVM = 0 
470 A$ = "a" 

480 INFUT "input target permanent dipole moments (y/n) " ; A$ 
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490 IF A$ = "Y" OR A$ = "y" THEN GOTO 510 

500 IF A$ = "N" OR A$ = "n" THEN GOTO 540 ELSE GOTO 480 

510 INRJT "permanent longitudineil moment in oersted-cmG (stem-to-bow +)"; PLM 
520 INHJT "permanent transverse mcament in oeirsted-onS (port-to-staxbcard +)"; 
PIM 

530 INPUT "permanent vertical rocinent in oersted-cm3 (downward +)"; PVM 
540 INPUT "target displacement in tons"; WT 
550 A$ = "a" 

560 INPUT "input target permeability coefficients or factors (c/f)"; A$ 

570 IF A$ = "C" OR A$ = "c" THEN GOTO 590 

580 IF A$ = "F" OR A$ = "f" IHEN GOTO 630 EISE GOTO 560 

590 INPUT "longitudinal permeability coefficient in cgs units"; KL 

600 INPUT "transverse permeability coefficient in cgs mits"; KT 

610 INPUT "vertical permeability coefficient in cgs units"; KV 

620 FL = KL / WT: FT = KT / WT: FV = KV / WT: GOTO 710 

630 INPUT "longitudinal displacement factor in cgs units"; FL 

640 INPUT "transverse displacement factor in cgs units"; FT 

650 INPUT "vertical displacement factor in cgs units"; FV 

660 KL = FL*WT: KT = FT*WT: KV = FV*WT 

670 A$ = "a" 

680 INPUT "generate a magnetic data file (y/n)"; A$ 

690 IF A$ = "Y" OR A$ = "y" THEN GOTO 710 

700 IF A$ = "N" OR A$ = "n" THEN GOTO 850 EI5E GOTO 680 

710 INPUT "magnetic data file name"; M$ 

720 ON ERROR GOTO 730: GOTO 740 

730 RESUME 710 

740 OPEN "0", #1, M$ 

750 WRITE #1, lAT, ING, DEC, DIP, HE, PIM, PIM, PVM, WT, KL, KT, KV, FL, FT, 

FV 

760 dOSE 
770 GOTO 850 

780 INPUT "magnetic data file name"; M$ 

790 ON ERROR GOTO 800: GOTO 810 

800 PRESUME 780 

810 OPEN "I", #1, M$ 

820 INPUT #1, lAT, ING, DEC, DIP, HE, PIM, PIM, PVM, WT, KL, KT, KV, FL, FT, 

FV 

830 CIDSE 

840 DECR = DEC * CON: DIPR = DIP * 0(34 
850 A$ = "a" 

860 INPUT "processing data entry by file/keyboard (f/k)"; A$ 

870 IF A$ = "F" OR A$ = "f" THEN GOTO 1450 

880 IF A$ = "K" OR A$ = "k" THEN GOTO 890 ELSE GOTO 860 

890 A$ = "a" 

900 INPUT "input saitpling period (y/n)"; A$ 

910 IF A$ = "Y" OR A$ = "y" THEN GOTO 930 

920 IF A$ = "N" OR A$ = "n" THEN GOTO 960 EISE GOTO 900 

930 INPUT "sairpling period in seconds"; DT 

940 IF DT <= 0 THEN PRINT : PRINT "must be greater than zero": PRINT : GOTO 
930 

950 GOTO 1090 
960 A$ = "a" 
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970 INHJT "ir^t maximum magnetic signal frequency (y/n)"; A$ 

980 IF A$ = "Y” OR A$ = "y" 1HEN GOTO 1000 

990 IF A$ = "N" OR A$ = "n" THEN GOTO 1030 EISE GOTO 970 

1000 INKIT "maximum magnetic signal frequency in Hertz"; M?0(F 

1010 IF MAXF <= 0 THEN PRINT : FKENT "must be greater than zero": PRINT : GOTO 

1000 

1020 GOTO 1080 

1030 INHJT '•minimum target slant range at CPA in meters"; MINRO 

1040 IF MINRO <= 0 THEN PRINT : KUOT "must be greater than zero": PEONT : 

GOTO 1030 

1050 INFUT "maximum magnetometer relative speed in knots"; MAXVMK 
1060 MAXVM = MAXVMK * KC»J: MAXF = 2 * MAXVM / MINRO 

1070 IF MAXVM <= 0 THEN IRINT : IRINT "must be greater than zero"; PRINT : 

GOTO 1050 

1080 DT = 1 / (2 * MAXF) : REM low pass filter, Nyquist sairpling rate 
1090 A$ = "a" 

1100 INRJT "irput integration time (y/n)"; A$ 

1110 IF A$ = "Y" OR A$ = "y" THEN GOTO 1130 

1120 IF A$ = "N" OR A$ = "n" THEN GOTO 1210 EISE GOTO 1100 

1130 INEUT "integration time in seconds"; IT 

1140 IF IT >= DT THEN GOTO 1170 

1150 PRINT ; HdlNT "IT = " + STR$(IT) + " seconds - minimum = " + STR$(DT) + " 
seconds": FRINT 
1160 GOTO 1130 

1170 NS = 2 * INT(IT / DT / 2) +1: REM adj number of samples per integration 
time 

1180 IF NS <= 151 THEN GOTO 1340 

1190 PRINT : FRINT "IT = " + STR$(IT) + " seconds - maximum = " + STR$(150 * 
DT) + " seconds": IRINT 
1200 GOTO 1130 

1210 INKJT "maximum target slant range at CPA in meters"; MAXRO 

1220 IF MAXRO <= 0 THEN ERINT : FRINT "must be greater than zero": FRINT : 

GOTO 1210 

1230 INFOT •’ininiinum magnetometer relative speed in knots"; MINVMK 
1240 MINVM = MINVMK * 

1250 IF MINVM <= 0 THEN PRINT ; FRINT "must be greater than zero": FRINT : 

GOTO 1230 

1260 IT = 2 * MAXRO / MINVM 
1270 IF IT >= DT THEN GOTO 1300 

1280 PRINT ; FRINT "IT = " + STR$(IT) + " seconds - minimum = " + STR$(DT) + " 
seconds": FRINT 
1290 GOTO 1100 

1300 NS = 2 * INr(IT / DT / 2) +1: REM adjusted number of samples per 

integration time 

1310 IF NS <= 151 THEN GOTO 1340 

1320 PRINT : PRINT "IT = " + STR$(IT) + " seconds - minintum = " + STR$(150 * 
DT) + " seconds": FRINT 
1330 GOTO 1100 
1340 A$ = "a" 

1350 INPUT "generate a processing data file (y/n)"; A$ 

1360 IF A$ = "Y" OR A$ = "y" THEN GOTO 1380 

1370 IF A$ = "N" OR A$ = "n" THEN GOTO 1510 ELSE GOTO 1350 
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1380 INRJT "processing data file name"; P$ 

1390 CM ERROR GOTO 1400: GOTO 1410 
1400 RESUME 1380 
1410 OPEN "O", #1, P$ 

1420 WRITE #1, DT, IT, NS 

1430 CIOSE 
1440 GOTO 1510 

1450 INPUT "processing data file name"; P$ 

1460 CM ERROR GOTO 1470: GOTO 1480 
1470 RESUME 1450 
1480 OPEN "I", #1, P$ 

1490 INPUT #1, DT, IT, NS 

1500 aosE 
1510 A$ = "a" 

1520 INPUT "kinematic data entry by file/keyboand (f/k)"; A$ 

1530 IF A$ = "F" OR A$ = "f" THEN GOTO 1790 

1540 IF A$ = "K" OR A$ = "k" THEN GOTO 1550 ELSE GOTO 1510 

1550 INPUT "magnetometer course in decimal degrees (0 if at rest)"; CM 

1560 INPUT "magnetometer speed in knots"; VMK 

1570 INPUT "magnetometer altitude in meters (below 0 is -)"; AM 

1580 INPUT "target course in decimal degrees"; CT 

1590 INPUT "target speed in knots"; VTK 

1600 INPUT "target depth in meters (above 0 is -)"; AT 

1610 CMR = CM * OON: CTR = CT * CX3N 

1620 WXK = VMK * SIN(CMR) - VTK * SIN(CTR) 

1630 WYK = VMK * COS (CMR) - VTK * OOS(CTR) 

1640 Z = AM + AT: REM vertical separation (- for magnetometer below target) 
1650 X = WXK: Y = WYK: GOSUB 4170 

1660 CR = T: C = CR / OON: REM magnetometer relative course 
1670 WK = R: REM magnetometer relative speed 
1680 A$ = "a" 

1690 INPUT "generate a kinematic data file (y/n)"; A$ 

1700 IF A$ = "Y" OR A$ = "y" THEN GOTO 1720 

1710 IF A$ = "N" OR A$ = "n" THEN GOTO 1860 EISE GOTO 1690 

1720 INPUT "kinematic data file name"; K$ 

1730 CM ERROR GOTO 1740: GOTO 1750 
1740 RESUME 1720 
1750 OPEN "O", #1, K$ 

1760 WRITE #1, CM, VMK, AM, CT, VTK, AT, Z, C, WK 
1770 dDSE 
1780 GOTO 1860 

1790 INPUT "kinematic data file name"; K$ 

1800 CM ERROR GCHO 1810: GOTO 1820 
1810 RESUME 1790 
1820 OPEN "I", #1, K$ 

1830 INPUT #1, CM, VMK, AM, CT, VTK, AT, Z, C, WK 
1840 CIDSE 

1850 CMR = CM * OCM: CTR = CT * CCM: CR = C * CON 

1860 MCMR = (CMR - DECR) : REM magnetometer magnetic course in radians 
1870 MCTR = (CTR - DECR) : REM target magnetic course in radians 
1880 HM = KL * HE * COS(DIPR) * COS(MCTR) : ITM = -KT * HE * COS(DIPR) * 
SIN(MCTR) 
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1890 IVM = KV * HE * SIN(DIPR) 

1900 CMX = (PIM + IIM) * SIN(MCrR) + (PIM + PIM) * aDS(MCTR) 

1910 EMY = (PIM + IIM) * OOS(MCTR) - (PIM + ITM) * SIN(MCrR) 

1920 EMV = PVM + IVM 
1930 X = CMX: Y = CMY 
1940 GOSUB 4170 

1950 CMIR = T: REM dipole azimuth relative to magnetic north 
1960 X = CM7; Y = R 
1970 GOSUB 4170 

1980 CM = R: CMR = T: REM dipole depression angle 
1990 CML = CMIR / 004: CM = CMR / 004 
2000 A$ = "a" 

2010 INHJT "generate a ccxtibined magnetic, processing & kinematic data file 
(y/n)"; A$ 

2020 IF A$ = "Y" OR A$ = "y" THEN GOTO 2040 

2030 IF A$ = "N" OR A$ = "n" IHEN GOTO 2210 ELSE GOTO 2010 

2040 U4HJT "ccxnbined magnetic, processing & kinematic data file name"; E$ 

2050 ON ERROR GCTTO 2060: GOTO 2070 
2060 RESUME 2040 
2070 OPEN "0", #1, E$ 

2080 WRITE #1, lAT, ING, DEC, DIP, HE, PIM, PIM, PVM, WT, KL, KT, KV, FL, FT, 
FV, DT, IT 

2090 WRITE #1, NS, CM, VMK, AM, CT, VTK, AT, Z, C, WK, IIM, TIM, IVM, CM, CML, 
OM, M$, P$, K$ 

2100 CLDSE 
2110 GCTO 2210 

2120 INHJr "combined magnetic, processing & kinematic data file name"; E$ 

2130 CM ERROR GOTO 2140: GOTO 2150 
2140 RESUME 2120 
2150 OPEN "I", #1, E$ 

2160 INHJT #1, lAT, I24G, DEC, DIP, HE, PIM, PIM, PVM, WT, KL, KT, KV, FL, FT, 
FV, DT, IT 

2170 INPUT #1, NS, CM, VMK, AM, CT, VTK, AT, Z, C, WK, IIM, ITM, IVM, CM, CML, 
CM, M$, P$, K$ 

2180 cmsE 

2190 DECR = DEC * CON: DIER = DIP * CON 

2200 CR = C * cm: CMIR = CML * Cm: CMR = CM * CCM 

2210 INPUT "required false alcomi rate in false alarms per hour"; FAR 

2220 PF = FAR * IT / 3600: REM false alarm probability 

2230 Y = PF: IF PF > .5 THEN Y = 1 - Y: REM inverse normal approximation 

2240 Y = SC2R(IOG(l / Y / Y) ) 

2250 Y = Y - (II + Y * (12 + 13 * Y)) / (1 + Y * (14 + Y * (15 + 16 * Y) ) ) 

2260 IF PF < .5 THEN Y = -Y 
2270 ZP = -Y 

2280 CHI = NS * (1 - 2 / 9 / NS + ZP * SQR(2 / 9 / NS) ) ‘ 3: REM inverse 

chi-square approximation 

2290 INPUT "magnetic noise in gamma"; SIG 

2300 INPUT "maximum lateral range in meters"; IRM 

2310 INPUT "lateral range step in meters"; ST 

2320 IF ST <= IRM THEN GOTO 2340 

2330 PRINT : PRINT "maximum step is " + STR$(IRM) + " meters": PRINT : GCTO 
2310 
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2340 NC = 2 * INT(LRM / ST) + 1: REM number of lateral range function values 
2350 IF NC <= 71 THEN GOTO 2370 

2360 EKENT : PRINT ''minirnum step is " + STR$(LRM / 35) + " meters": PRINT : 
GOTO 2310 

2370 ATT = DT * NS: REM adjusted integration time 

2380 W = WK * KCW: REM magnetometer relative speed in meters/second 

2390 DS = W * DT: REM distance between samples on the relative track 

2400 XO = -(NC - 1) / 2 * ST 

2410 FOR I = 0 TO NC - 1 

2420 X0(I) = XO 

2430 X = XO: Y = Z 

2440 GOSUB 4170 

2450 RO = R: R0(I) = R: REM target slant range at CPA in meters 

2460 DELR = T: REM target depression angle complement at CPA in radians 

2470 IF RO = 0 THEN GOTO 2800: REM zero lateral range and vertical separation 

2480 CMF = EM / 10 / RO ‘ 3: REM dipole mcment factor 

2490 OylR = CR - DECR: REM target relative magnetic course 

2500 L = OOS(CMR) * COS(CMR - OMLR) 

2510 M = COS (DELR) * OOS(CMR) * SIN(CMR - CMIR) - SIN(DEER) * SIN(CMR) 

2520 N = -SIN(DEnJR) * COS(CMR) * SIN(CMR - CMIR) - COS(DEIR) * SIN(CMR) 

2530 LI = OOS(DIPR) * OOS(CMR) 

2540 Ml = COS(DEIR) * OOS(DIPR) * SIN(CMR) - SIN(DEIR) * SIN(DIPR) 

2550 N1 = -SIN(DEIR) * OOS(DIPR) * SIN(CMR) - COS (DEER) * SIN(DIRR) 

2560 A2 = 2*L*L1-M*M1-N*N1: REM Anderson Function Coefficient 
2570 A1 = 3 * (N * LI + L * Nl) : REM Anderson Function Coefficient 
2580 A0 = 2*N*N1-L*L1-M* Ml: REM Anderson Function Coefficient 
2590 SUM = 0: HMAX = 0: HMIN = 0 
2600 FOR J = 0 TO NS - 1 

2610 S = (J - (NS - 1) / 2) * DS: BA = S / RO: E®I Anderson Function Argument 

2620 AF = 1/ (1 + BA* BA) ' 2.5: REM Anderson Function Factor 

2630 HSF = (A2 *BA*BA + A1*BA + AO) * AF: REM magnetic signal factor 

2640 HS(I, J) = EMF * HSF: REM magnetic signal value 

2650 IF HS(I, J) > HMAX THEN HMAX = HS(I, J) 

2660 IF HS(I, J) < HMEN THEN HMEN = HS(I, J) 

2670 SUM = SUM + HSF * HSF 
2680 NEXT J 

2690 H(I) = HMAX - HMEN 
2700 K(I) = SQR(SUM) 

2710 W = -ZP + EMF * SQR(SUM) / SIG 

2720 lAM = EMF * EMF * SUM / (SIG * SIG) : AA = NS + lAM: BB = 1 + lAM / (NS + 
lAM) 

2730 ZN = -SQR(2 * CHI / BB) + SQR(2 *AA/BB-1): X1=W 

2740 GOSUB 4230 

2750 IF Y1 > 1 THEN Y1 = 1 

2760 HXX:(I) = Yl: XI = ZN 

2770 GOSUB 4230 

2780 IF Yl > 1 THEN Yl = 1 

2790 roSL(I) = Yl 

2800 XO = XO + ST 

2810 NEXT I 

2820 C = C / 360: C = (C - INT(C)) * 360 
2830 IF C < 0 THEN C = 360 + C 
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2840 OIL = (CML + DEC) / 360: CML = (CML - INTCCML)) * 360 
2850 IF CML < 0 THEN CML = 360 + 

2860 PRINT : A$ = "a” 

2870 INPUT "generate a program data file (y/n)"; A$ 

2880 IF A$ = "Y" OR A$ = "y" THEN GOTO 2900 

2890 IF A$ = "N" OR A$ = "n" THEN GOTO 3230 ELSE GOTO 2870 

2900 INPUT "program data file name"; D$ 

2910 C»J ERROR GOTO 2920: GOTO 2930 
2920 RESUME 2900 
2930 OPEN "O", #1, D$ 

2940 WRITE #1, lAT, WG, DEC, DIP, HE, PIM, PIM, PVM, WT, KL, KT, KV, FL, FT, 
FV, DT, IT, ATT 

2950 WRITE #1, NS, CM, VMK, AM, CT, VTK, AT, IIM, TIM, IVM, Z, C, WK, EM, CML, 
CM, FAR, PF, SIG, ST 

2960 WRITE #1, LRM, DS, NC, M$, P$, K$, E$ 

2970 FOR I = 0 TO NC - 1 

2980 WRITE #1, X0(I), PDCC(I) , PDSL(I) , K(I) , H(I) , R0(I) 

2990 NEXT I 

3000 FOR I = 0 TO NC - 1 
3010 FOR J = 0 TO NS - 1 
3020 WRITE #1, HS(I, J) 

3030 NEXT J 
3040 NEXT I 
3050 CLOSE 
3060 GCTO 3230 

3070 INPUT "program data file name"; D$ 

3080 CM ERROR GOTO 3090: GOTO 3100 
3090 RESUME 3070 
3100 OPEN "I", #1, D$ 

3110 INPUT #1, lAT, ING, DEC, DIP, HE, PIM, PIM, PVM, WT, KL, KT, KV, FL, FT, 
FV, DT, IT, ATT 

3120 INPUT #1, NS, CM, VMK, AM, CT, VTK, AT, IIM, TIM, IVM, Z, C, WK, EM, CML, 
CM, FAR, PF, SIG, ST 

3130 INPUT #1, IRM, DS, NC, M$, P$, K$, E$ 

3140 FOR I = 0 TO NC - 1 

3150 INPUT #1, X0(I), roCC(I), PDSL(I) , K(I) , H(I) , R0(I) 

3160 NEXT I 

3170 FOR I = 0 TO NC - 1 

3180 FOR J = 0 TO NS - 1 

3190 INPUT #1, HS(I, J) 

3200 NEXT J 
3210 NEXT I 
3220 dOSE 

3230 PRINT : A$ = "a" 

3240 INPUT "print encounter parameter values (y/n)"; A$ 

3250 IF A$ = "Y" OR A$ = "y" THEN GOTO 3270 

3260 IF A$ = "N" OR A$ = "n" THEN GOTO 3780 ELSE GOTO 3240 



3270 LPRINT 

3280 LFRINT "program file name " + N$ 
3290 LFRINT "program data file name " + E»$ 
3300 LFRINT "magnetic data file name " + M$ 
3310 LFRINT "processing data file name " + P$ 
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3320 LFRINT 
3330 LEKDW 
3340 IfKENT 
lAT 

3350 LPRINT 
IWG 

3360 LPEmrr 
DEC 

3370 LPRINT 
DIP 

3380 LEKCNT 
HE 

3390 LPRnn’ 
PLM 

3400 LPRINT 
PIM 

3410 LPRINT 
PVM 

3420 LPRINT 
WT 

3430 LPRINT 
KL 

3440 LPRINT 
KT 

3450 LPRINT 
KV 

3460 LPRINT 
FL 

3470 LPRINT 
FT 

3480 LEKENT 
FV 

3490 LPRINT 
DT 

3500 LPRINT 
IT 

3510 LPRINT 
AIT 

3520 LEKENT 
NS 

3530 LEKENT 
CM 

3540 LPRINT 
VMK 

3550 LEKENT 
AM 

3560 LEKENT 
CT 

3570 LPRINT 
VTK 

3580 LPRINT 
AT 



'kinematic cJata file name 

'cxjmbined magnetic, processing & kinematic data file name 
'encounter latitude (decimal degrees) " 

'encounter longitude (decimal degrees) " 

'encounter variation (decimal degrees) " 

•encounter dip cingle (decimal degrees) " 

'encounter magnetic field intensity (gamma) " 

'permanent longitudinal moment (oersted-arG) " 

'permanent transverse mcroent (oersted-crti3) " 

'permanent vertical moment (oersted-cm3) " 

'target displacement (tons) " 

'target longitudinal permeability coefficient " 

•target transverse permeability coefficient " 

•target vertical permeability coefficient " 

•target longitudinal permeability factor •• 

•target transverse permeability factor •• 

•target vertical permeability factor •• 

•sampling period (seconds) •• 

•integration time (seconds) •• 

•adjusted integration time (seconds) •• 

•number of samples per encounter •• 

•magnetometer course (decimal degrees) •• 

•magnetometer speed (knots) •• 

•magnetometer altitude (meters) •• 

•target course (decimal degrees) •• 

•target speed (knots) •• 

•target depth (meters) •• 



" + K$ 
•• + E$ 
SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 

SPC(2) 
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3590 LFRINT "magnetiOneter relative course (decimal degrees) " 

C 

3600 LFKENT "magnetometer relative ^>eed (knots) " 

WK 

3610 LPEUNT "magnetcmeter-target vertical s^jaration (meters) " 

Z 

3620 LFKENT "target induced longitudinal dipole moment (oersted-cm3) " 
IIM 

3630 LFKENT "target induced transverse dipole moment (oersted-cm3) " 
TIM 

3640 LFKENT "target induced vertical dipole moment (oersted-cm3) " 
IVM 

3650 LFKENT "magnetic dipole mcment (oersted-cm3) " 

EM 

3660 LFKENT "dipole moment azimuth (decimal degrees) " 

CML 

3670 LFKENT "dipole moment depression angle (decimail degrees) " 

CM 

3680 LFKENT "distance between samples on the relative track (meters)" 
DS 

3690 LFKENT "false alarm rate (false alarms per hour) " 

FAR 

3700 LFKENT "false alarm probability " 

PF 

3710 LFKENT '•magnetic noise (gamma) " 

SIG 

3720 LFKENT "maximum lateral range (meters) " 

LRM 

3730 LFKENT "lateral range step (meters) " 

ST 

3740 LFKENT "nuniber of lateral range function values " 

NC 

3750 FOR I = 0 TO 15 
3760 LFKENT 
3770 NEXT I 
3780 FKENT : A$ = "a" 

3790 INFUT "print lateral range function values (y/n) " ; A$ 

3800 IF A$ = "Y" OR A$ = "y" THEN GOTO 3820 

3810 IF A$ = "N" OR A$ = "n" THEN (SOTO 3900 ELSE GOTO 3790 

3820 LFKENT D$; " lateral range function values" 

3830 LFKENT : LFKENT 

3840 LFKENT "L p(o^) P(sl) H R 

k" 



SPC(2) 

SFC(2) 

SPC(2) 

SF>C(2) 

SPC(2) 

SPC(2) 

SFC(2) 

SPC(2) 

SF>C(2) 

SFG(2) 

SPC(2) 

SPC(2) 

SE>C(2) 

SF>C(2) 

SF>C(2) 

SFC(2) 



3850 LFKENT "meters gamma meters" 

3860 LFKENT 

3870 FOR I = 0 TO NC - 1 

3880 LFKENT X0(I); TAB(IO) ; PDCC(I) ; TAB(24) ; PDSL(I) ; TAB(38) ; H(I) ; TAB(52) ; 
R0(I); TAB(70); K(I) 

3890 NEXT I 



3900 FKENT : A$ = "a" 

3910 INRJT "print magnetic signal values (y/n)"; A$ 
3920 IF A$ = "Y" OR A$ = "y" THEN GOTO 3940 
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3930 IF A$ = "N" OR A$ = "n" THEN GOTO 4120 EISE GOTO 3910 

3940 FKENT : PRINT "identify signal by encounter lateral range index" 

3950 PRINT "lateral range equals the index times " + STR$(ST) + " meters" 

3960 PRINT "index values: + STR$((NC - 1) / 2) + " to " + STR$((NC - 1) / 

2) 

3970 INPUT "lateral range index"; K: I = K + (NC - 1) / 2 

3980 LPRINT : LFRINT : LPRINT 

3990 LERINT D$; " magnetic signal values" 

4000 LPRINT "for a lateral range of "; K * ST; " meters" 

4010 LERINT 

4020 LERINT "relative CPA distance"; TAB(35) ; "magnetic signal" 

4030 LERINT "s in meters"; TAB(35) ; "Hs in gamma" 

4040 LERINT 

4050 FOR J = 0 TO NS - 1 

4060 LERINT (J - (NS - 1) / 2) * DS; TAB(35) ; HS(I, J) 

4070 NEXT J 

4080 PRINT : A$ = "a" 

4090 INIUT "print magnetic signal values for a different lateral range (y/n)" 
A$ 

4100 IF A$ = "Y" OR A$ = "y" THEN GOTO 3940 

4110 IF A$ = "N" OR A$ = "n" THEN GOTO 4120 EIEE GOTO 4090 

4120 PRINT : A$ = "a" 

4130 INPUT "continue to use the program (y/n)"; A$ 

4140 IF A$ = "Y" OR A$ = "y" THEN GOTO 10 

4150 IF A$ = "N" OR A$ = "n" THEN GOTO 4160 EIEE GOTO 4130 

4160 END 

4170 R = SQR(X * X + Y * Y) : REM rectangular to polar conversion 
4180 IF R = 0 THEN T = 0: RETORN 

4190 IF ABS(X / R) = 1 THEN Q = SGN(X) * (PI / 2) ELSE Q = ATN(X / R / SQR(1 
X * X / R / R) ) 

4200 IF ABS(Y / R) = 1 THEN T = (PI / 2) * (1 - SOI(Y)) ELSE T = (PI / 2) - 
ATN(Y / R / SQR(1 - Y*Y/R/R)) 

4210 IF Q < 0 THEN T = 2 * PI - T 
4220 RETURN 

4230 Y1 = XI: IF XI < 0 THEN Y1 = -Yl: REM normal approximation 
4240 G = 1 / (1 + QO * Yl) 

4250 Yl = EXP(-Y1 * Yl / 2) / SQR(2 * PI) * G * (Q1 + G * (Q2 + G * (Q3 + G * 
(Q4 + G * Q5) ) ) ) 

4260 IF XI >= 0 THEN Yl = 1 - Yl 
4270 RETURN 
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